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Abstract 

We develop a numerical approach based on our recent analytical model of fiber structure in the left ventricle of the human 
heart. A special curvilinear coordinate system is proposed to analytically include realistic ventricular shape and myofiber 
directions. With this anatomical model, electrophysiological simulations can be performed on a rectangular coordinate grid. 
We apply our method to study the effect of fiber rotation and electrical anisotropy of cardiac tissue (i.e., the ratio of the 
conductivity coefficients along and across the myocardial fibers) on wave propagation using the ten Tusscher-Panfilov 
(2006) ionic model for human ventricular cells. We show that fiber rotation increases the speed of cardiac activation and 
attenuates the effects of anisotropy. Our results show that the fiber rotation in the heart is an important factor underlying 
cardiac excitation. We also study scroll wave dynamics in our model and show the drift of a scroll wave filament whose 
velocity depends non-monotonically on the fiber rotation angle; the period of scroll wave rotation decreases with an 
increase of the fiber rotation angle; an increase in anisotropy may cause the breakup of a scroll wave, similar to the mother 
rotor mechanism of ventricular fibrillation. 
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Introduction 

The modeling of cardiac electrical function is a well-established 
area of research that began with early models of cardiac cells 
developed by D. Noble [1]. 

The importance of modeling in cardiology comes from the 
widespread prevalence of cardiac disease. For example, sudden 
cardiac death is the leading cause of death in the industrialized 
world, accounting for more than 300,000 victims annually in the 
US alone [2]. In most cases, sudden cardiac death is a result of 
cardiac arrhythmias that occur in the ventricles of the human 
heart [2]. 

When studying cardiac arrhythmias, it is important to 
understand that they often occur at the level of the whole organ 
and in these situations cannot be reproduced in single cells. 
Therefore, it is very important to model cardiac arrhythmias at the 
tissue level, preferably using an anatomically accurate represen- 
tation of the heart. Compared to modeling at the single-cell level, 
anatomical modeling started much more recently [3,4]. Using 
anatomical models, researchers have been able to obtain 
important results on the 3-D organization of cardiac arrhythmias 



in animal [5] and human [6] hearts. Moreover, the defibrillation 
process has been investigated [7], and the effects of mechano- 
electrical coupling on cardiac propagation have recently been 
modeled [8, 9]. Multi-scale anatomical cardiac modeling is becom- 
ing increasingly prominent in medical and pharmaceutical 
research [10]. 

In a broad sense, an anatomical model of the heart is a 
combination of models of cardiac cells and anatomical data. The 
development of models of the electrical and mechanical functions 
of cardiac cells is a well-established area of research, and many 
models have been developed, including models of the human 
cardiac cells [11-20]. The anatomical data necessary for cardiac 
modeling include not only the heart's geometry but also its 
anisotropic properties. Although the geometry of the heart can be 
obtained from routine clinical procedures such as MRI or CT 
scans [21,22], anisotropy data are much more challenging to 
acquire. Currently, the acquisition can be done on explanted 
hearts only, using either direct histological measurements or time- 
demanding DT-MRI scans [23-26]. In addition to experimental 
noise, even perfect measurements will grant only the particular 
anisotropy of the imaged heart. Thus, to study the effects of 
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anisotropy on wave propagation, one needs to vary the anisotropic 
properties and to separate the anisotropy effects from other 
factors. All these questions can be addressed with the development 
of models that account for the anisotropy of the heart using 
analytical or numerical tools. 

In a previous article [27], we described an axisymmetric model 
of the left ventricle (LV) of the human heart. In the model, we 
represented the LV shape (including positions of cardiac fibers) as 
analytical functions of special curvilinear coordinates defined on a 
rectangular domain. Our model allowed the generation of not only 
a default architecture of anisotropy closest to the reality but also 
intermediate architectures that can be used to study the effects of 
any specific element of anisotropy on wave propagation in the 
heart. 

In this paper, we build on our previous approach in two ways. 
First, we develop a numeral scheme for the integration of 
equations for wave propagation in our anatomical model of the 
LV, which is the best possible way to account for anisotropy. In 
particular, we develop a model on a rectangular domain and 
represent anisotropy and the LV shape by means of parameter 
changes. Second, we vary the geometry and anisotropy parameters 
to study how the rotation of the fiber orientation affects wave 
propagation and show that rotational anisotropy accelerates the 
spread of electrical excitation in the heart. We also study the 
behavior of scroll waves and their filaments. We show that the 
scroll waves drift and we calculate their drift velocity and period of 
rotation depending on the fiber rotation angle and the diffusion 
coefficients ratio. 

Model Description 

Geometrical model of the LV 

In our model, the LV is represented as a body of revolution 
around the vertical axis Oz with the shape fitted to experimental 
data (for details, see [27]). A section of the LV is shown in Fig. 1. 
The rotation of the blue line delineates the epicardial surface, 
while the rotation of the red line yields the endocardial surface of 
the LV. 

In our model, each point of the LV has three local coordinates 
(y, ij/, (j)). The coordinate y (yo^y^yi) gives us points between the 
endo- and epicardial surfaces in Fig. 1, i.e., it is a measure for 
transmural depth; the coordinate xj/ (0<i/^<7r/2) is explained in 
Fig. 1; and the coordinate (j) (0<(^<27r) is the rotation angle 
around the vertical axis Oz- The local coordinates are linked with 
the cylindrical coordinate system (CS) (p, ^, z) by the following 
formulae: 

P = {n + {\-y)l)eM), (1) 
z = (ZA + ( 1 - y)h) ( 1 - sin ^) + (y - 1 )h, (2) 

with: 

ec = ecosi/^ + (l— e)(l— sini/^). 
Below, we also use 

e^ = esini/^ + (l — e) cosi/^. 

The physical meaning of the parameters is as follows: fb is the LV 
cavity radius on the LV equator, Zh is the LV cavity depth, / is the 




Figure 1. A radial section of the endocardial (solid red line) and 
epicardial (dashed blue line) surfaces of the LV model, from 

[27]. 

doi:10.1371/joumal.pone.0093617.g001 

LV wall thickness on the LV equator, h is the LV wall thickness on 
the LV apex, and ee [0,1] is a dimensionless parameter influencing 
the conicity-ellipticity characteristic of the LV shape. In this 
system, y — yo gives the epicardium and y — y\ gives the endocar- 
dium; the value i/^ = 0 describes the upper (basal) boundary of the 
LV. 

Following Pettigrew's idea [28] about spiral surfaces and 
semicircle chords mapping on the surfaces, our LV model consists 
of spiral surfaces on which a set of curves is defined. The details 
are described in our previous work [27] and are summarized in 
Appendices A (spiral surfaces construction) and B (fiber equations). 
In Figs. 2 and 3, we show common views of spiral surfaces and the 
fibers inside them. 

We demonstrated in [27] that the model approximates the real 
fiber orientation field in the LV reasonably well. A comparison of 
true fiber angle a and helix angle ai with MRI data showed that 
fiber architecture in the equatorial region of the heart was well 
reproduced in our anatomical model. In the middle (by height) 
and apical areas, the angles were reproduced both qualitatively 
and quantitatively well; the difference between the model and the 
experimental data was not more than 25°. 

Overall, we can consider the anatomical model as a map from a 
rectangular domain yo— ? — 0<i/^<7r/2; 0<(^<27r to the shape 
of LV with anisotropy explicitly given by Eqs. (B.1)-(B.3). The 
total fiber rotation angle Aai is defined as the difference between 
the epicardial and endocardial helix angles ai measured at the LV 
basal zone i/^ = 7r/8 (see [29] for details). It can be varied by 
changing the values of the parameters yo ^nd yi'. 
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o^i(y*,'A*) = (P$rnv)|y=y*,^=^*,0=o, 
AaiCyoJi) = ai(yi,7i/8)- ai(yo,7r/8), 

where p is the tangent vector to the parallel y = const, ij/ = const at 
the point y = y , ij/ = ij/ , (j) = 0 (the geometric model is axisymmet- 
ric; therefore the choice of 0 is arbitrary); n is the normal vector to 
the epicardium passing through the same point; v is the fiber 
direction vector (it is defined by Eq. (B.1)-(B.3)); pr is the 
projection operator; and (u,w) denotes the angle between the 
vectors u and w. We have chosen the value i/^ = 7r/8 because in 
this case, the total fiber rotation angle changes uniformly enough 
depending on the values of yo,i we use (see Table 1 and section 
"Parameter values" below). For short, below we will denote 
Aai(yo5 yi) as a without argument. We use it in the present study to 
investigate the effect of fiber rotation on wave propagation in the 
heart. 

Electrophysiological model 

To describe the excitation of cardiac tissue, we use the detailed 
ionic model for human ventricular cells from [6,11]. The model 
uses reaction-diffusion equations to describe the evolution of the 
transmembrane potential u = u{r, t): 

J=div(Dgrad«)-^, (4) 



+ Iks +lKl+Ito+ Ino + IbNa + 1 Cat + . , 

(5) 

hCa + InqK + InuCu + IpCa + IpK • 

Here, the intracellular processes are captured by li^n — lion^^ ^ 
which is the sum of the ionic transmembrane currents; is the 
capacitance of the cell membrane. The locally varying diffusion 
matrix D accounts for myofiber anisotropy. As in [6] , the diffusion 
matrix D = (D^ was computed using the following formula 

D''=D25ij + (D^-D2)viVj, (6) 




Figure 2. A spiral surface. The lines on the surface have equations 
p = const and = const. Color corresponds to height (z coordinate). 
doi:10.1371/journal.pone.0093617.g002 




Figure 3. A spiral surface viewed from the top (left panel) and 
side (right panel). Two myofibers are displayed in red and black. 
doi:10.1371/journal.pone.0093617.g003 

where Di and D2 are the diffusion coefficients along and across the 
fibers, V is the unit vector of fiber direction, j are Cartesian 
indices, and djj is the Kronecker symbol. 

Numerical integration scheme and boundary conditions 

The aim in this paper is to present a numerical procedure that 
allows us to use the analytical representation of cardiac anatomy 
and anisotropy described in the previous section. In particular, as 
our anatomical model is just a map from a rectangular domain 
yQ<y<yj; 0<i/^<7r/2; 0^(j)<2n to the shape of LV, we can 
formulate our approach in that rectangular domain. The shape of 
the heart, as well as anisotropy, will then be a curvilinear 
coordinate system (l)-(2) defined on that domain. We need to 
recalculate Eq. (4) with no-flux boundary conditions in those 
coordinates. A long computation presented in appendix C results 
in the following expression of the diffusion term: 

d.v(D grad «) = _ + _ , (7) 

where = y, ^2 - 'A? ^3 - and/>yt and q^i are coefficients given by 
the explicit analytical Eqs. (G.12) and (C.13) that depend only on 
the geometry of the LV and on the diffusion matrix. This 
representation is similar to that presented by Sridhar et al. in [30] . 
However, in our work, it was done for the 3-D case and for a 
special form of the diffusion matrix (see Eq. (6)). 

For numerical integration of the model (B. 1)-(B.3), (4), (5), we 
use the explicit finite difference method on a discrete grid in the (y, 
(j)) space. We initially use a uniform grid with the y-indices of 
the nodes denoted as z = 0, 1, ... Uy', the i/^-indices asj = 0, \, ... n^; 
and the (/>-indices as /:=0, 1, ... 

Although this grid is uniform in the (y, ij/, (j)) space, it is non- 
uniform in real space because distances between the grid points 
substantially decrease when ij/ approaches n/2, which is similar to 
the situation at the pole in a polar CS. Note, however, that even 
the cubic Cartesian lattice in real space is not uniform with respect 
to anisotropic diffusion due to the curved space interpretation of 
anisotropy [31,32]. To account for this problem, we exclude some 
points from our uniform grid in the following way. We first choose 
a threshold value of distance 4nin- Then, at y — yi (i.e., at the 
epicardial surface) and any given ij/ = ij/j, we calculate the distances 
between the node at (/> = 0: (y = yi, lA = 'Aj? ^-^) node (y = yi, 
^ — (j) — (l)k). We find minimal k satisfying two conditions: (1) the 
distance to the ^-node from the node at (/> = 0 is more than the 
threshold value ^/minj and (2) is a divisor of n^. We denote this 
number as Kj (as it depends on xj/j). Ifij/ is far from n/2, then Kj— 1 
and we use all nodes of our uniform (y^, ij/j^ grid. When ^ 
approaches the value of 7i/2, Kp> 1 and we drop all nodes between 
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Table 1. 


Dependence of the total fiber angle on the model parameters yo,i- 
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The helix angle near the base at 
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the epicardium the endocardium 


rotation angle a 
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0.3 


0.55 


-13° 3° 


16° 
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0.2 


0.7 


-40° 29° 


69° 


3 


0.1 


0.85 


-69° 64° 


133° 
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-87° 87° 


174° 
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(y -yi,^- ^j, <^ = 0) and (y = yi, lA = 'Ajj <^ = ^5 then the next node 
will be (y = yi, ^-'^p ^ = 2Kj) etc. Thus, only nodes with the ^- 
indices 0, AJ, . . . will be taken for evaluation of the Laplacian. 
After each time step, we compute values for all variables of our 
model in the omitted nodes using linear interpolation. With this 
approach, we reduce the number of grid elements in the (j)- 
direction in the apical region; otherwise, we would have needed to 
lower the maximal time step in our explicit integration scheme, 
which would have slowed our computations significantly. 
The no-flux boundary condition in our problem is 



nD grad w = 0, 



(8) 



where n is a normal to surface. We rewrite this equation in our 
special CS and obtain the following expression: 



du du du ^ 



(9) 



where Cy, c^j^, C(j, are coefficients given by Eqs. (D.14) and (D.16) in 
Appendix D. In order to satisfy the boundary condition, we add 
nodes at the domain boundaries in the following way. In the 
special GS, the domain of integration is a rectangle (with periodic 
boundary in (/>), and at 1^ = 71/2 we have a pole (i.e., also no 
boundary). Therefore, we have only three boundaries, namely at 
y = yo (i-e., ^ = 0), y = yi [i-n^, and 1^ = 0 (/'=0). Fictitious layers 
with (— 1, J, A:), [uy +1, J, ^), and (z, — 1, A:) are added. We then solve 
equation (9) on the three boundary surfaces to find values in the 
added nodes. Subsequently, we can compute Laplacian at all other 
nodes in the domain using, if necessary, the values in the 
additional nodes. Due to this procedure, all the nodes lying on the 
boundaries will satisfy the boundary conditions. 

We have programmed this approach using C in the GodeBlocks 
IDE, a Mingw compiler. The calculations were performed in 
operating systems Windows 7 and Linux. The OpenMP and MPI 
technologies have been used for parallelization, and Paraview and 
Irfan have been used for visualisation. The formal parameters of 
our numerical scheme have been given in the methods section 
above. Such an approach allows us to compute various regimes of 
wave propagation in a model of LV with good representation of 
boundary conditions and to study various effects of anisotropy on 
wave propagation patterns. 

We also studied the dynamics of scroll waves using the ten 
Tusscher-Noble-Noble-Panfilov (TNNP) model [11] and various 
anisotropy parameters. We initiated a scroll wave using the S1-S2 
stimulation protocol and studied its dynamics 12 s. For drift 
velocity and rotation period calculations, we take into account only 
the last 8 seconds of the simulations to exclude transient processes. 
We calculated average periods in the section = 0, that is, x — p> 



0, J = 0. Filaments were analyzed as reported in [6,33]. Finally, we 
computed their average velocity v in the Cartesian CS. 

Parameter values 

We used the following parameters from Table 2 in [29]: the LV 
equatorial radius = 23 mm, the equatorial wall thickness 
f = 10 mm, the LV cavity depth z^ = 53 mm, the apical wall 
thickness h^=7 mm, the conicity-ellipticity parameter e = 0.9, and 
the spiral surface torsion angle (^^^^ = Stt (see Fig. 3c in [29]). The 
threshold distance between the adjacent nodes a'min was set to 
0.3 mm. Currently, the first four parameters are measured using 
modern experimental techniques such as MRI (see [34-36]). The 
values used in our paper are in agreement with these experimental 
data. 

Our mesh had a distance of 0.2 to 0.3 mm between the nodes 
and before the deletion of the nodes described above had Uy = 40, 

= 300, — 800. The diffusion coefficient along the fibers was 
Z)i=0.3 mm^/ms. The diffusion coefficient across the fibers D2 
was varied between different experiments depending on whether 
we modeled isotropy or anisotropy. 

We applied our approach to study the effect of anisotropy on 
the spread of excitation in the heart. In particular, we initiated a 
wave at several locations and studied how the wave arrival time 
depends on the two main features of anisotropy. Our first 
anisotropy parameter was the ratio of the diffusion coefficients 
Di'.D2 along and across the fibers. Also, we independently varied 
the total rotation angle of fibers through the myocardial wall by 
adjusting yo? ?! ^^d keeping wall thickness constant. We also 
compared our results with the spread of excitation in an isotropic 
model of the LV where Di — D2. 

The dependency of wave velocity on the direction of 
propagation in the heart was measured in [37-39]. Experimental 
data show that the ratio of longitudinal to transverse conduction 
velocities ranges between 3 [38] and 2.1 [39]. Since D cx c^, where 
c is wave velocity (see, e.g., [40]), we used the ratios 
Z)i:Z)2 = 1:0.1 1 1 and Z)i:A = 1:0-25, which correspond to the 
experimental data. These anisotropy ratios were also used in the 
modeling studies [33,41-44]. 

For point stimulation, we increased the value of the variable u 
from the resting potential of —86.2 mV to u — 0 mV at the first 
time step in small regions located at three different sites. In the A 
series, it was a small epicardial region at the apex; in the B series, 
at the centre of LV epicardium; in the C series, at the centre LV 
endocardium (see Table 2). 

In this paper, we study the effect of the fiber rotation on the 
spread of excitation. With this purpose, we generated a series of 
LV models that differ in the total fiber rotation angle through the 
myocardial wall. The parameters of the model are listed in Table 1. 
Note that although the values of yo and yi differ between the 
models, they affect only fiber rotation, and the LV geometry is 
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Table 2. The initial excitation areas. 
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i>ny-4 


^•-(n^/2)|<2 


/c<4 


central epicardium 
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^•-(n^/2)|<2 


/c<4 


central endocardium 
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exactly the same for all models due to the rescaling procedure 
described in Sec. 1 .4 in our previous article [27] . Also note that the 
change in fiber rotation results in the change in fiber angle at the 
epicardial surface (see column "a" of Table 1). 

The time step was equal to 0.005 ms for the isotropic cases and 
0.01 ms for the anisotropic ones. 

We considered that a wave came to a node when the potential 
in the node was more than —80 mV the first time. 

Numerical Results 

Activation nnaps 

Figure 4 shows the wave arrival time after stimulation of the 
small apical epicardial zone A for ratios Di\D2 equal to 1:0.1 1 1 or 
1:0.25 (shown at the top of the figure) and four different rotation 
angles of the fibers, which are displayed in the left column. We see 
that in this first example, all the figures are axisymmetrical, which 
is a consequence of the axisymmetric properties of our model and 
the initial conditions. 

We observe that for the low rotation angle (the upper row), the 
speed of the upward wave propagation for the diffusion coefficient 
ratio of 1:0.111 is substantially smaller than that for the ratio of 
1:0.25. However, if the fiber rotation angle increases (the lower 
rows), the difference in the speed between the two anisotropy 
ratios decreases. For the rotation angle of 1 74° (the lower row), the 
excitation patterns for both anisotropy ratios become close to each 
other. Thus, we observe that fiber rotation increases the velocity of 
the spread of excitation and also decreases the effect of anisotropy. 

Let us now consider the case of lateral stimulation for a given 
ratio of the diffusion coefficients of 1:0. 1 1 1; the results for epi- and 
endocardial stimulation are presented in Figs. 5 and 6. After 
epicardial stimulation, the wave initially follows the fiber direction. 
In Fig. 5, we note a displacement of the early activation zones (red) 
due to the change in fiber direction at the epicardium in our 
anatomical model (see column 4 of Table 1). However, for 
endocardial stimulation (Fig. 6, the first column), the shift of the 
early activation zone (red) on the epicardium is attenuated by fiber 
rotation. As in Fig. 5, we see that an increase in the rotation angle 
causes a decrease in the arrival time. In addition, in the second 
column in Fig. 5, the excitation patterns have a clear V shape on 
the surface opposite the stimulation site, which is a mere 
consequence of the shape of the heart (see [3,4]). 

Figs. 7 and 8 show the arrival time after lateral stimulation in a 
case of decreased anisotropy, i.e., for a diffusion coefficient ratio of 
1:0.25. The excitation patterns resemble those from Fig. 5 and 
Fig. 6 respectively, with similar effects of fiber rotation on the 
epicardial stimulation and V-shaped patterns. Here, we also 
observe that an increase in the rotation angle decreases the overall 
excitation time. A compensating effect of fiber rotation on the 
degree of anisotropy can be noted. The difference between the 
corresponding panels in Fig. 5 and Fig. 7 (and also Fig. 6 and 
Fig. 8) is more pronounced for a lower fiber rotation rate (rotation 
angle 16°). 



Average speed of excitation 

Now let us quantify the effects of rotation and anisotropy on 
wave propagation. In order to do this, we use the following 
procedure. We group all points of the heart to bins differing by 
their "distance" from the stimulation point. We define the distance 
as the arrival time from the stimulation to a given point. To 
calculate the distances, we perform simulations in which we 
initiate a wave at the same locations as in Figs. 4-8. However, for 
the isotropic medium, we use a diffusion coefficient of 
Z)i=Z)2 = 0.3 mm^/ms. We generate 40 groups in which points 
differ in the arrival time by 2 ms. Then we determine the average 
arrival time for each of these groups for various anisotropic 
conditions and compare these arrival times to the arrival time in 
the isotropic model. As in the isotropic model, the velocity of wave 
propagation in all directions is the same; this dependence gives us 
the dependence of the wave arrival time on the distance from the 
stimulation point. 

In Fig. 9, the red lines correspond to a = 1 74° and the black 
lines correspond to a = 16° fiber rotation angle in the LV wall. The 



total fiber 

rotation Di:D.^l:0.111 
allele a 



Di:D2=l:0.25 



16° 




69° 



133° 



1163 




174° 



Figure 4. Arrival times, in ms, of the waves after point 
stimulation at the apex for various values of anisotropy and 
fiber rotation. The values of anisotropy are shown at the top of the 

figure and the values of the fiber rotation are shown in the left column. 
For details, see Table 1. Arrival times are color-coded in ms. 
doi:1 0.1 371/journal.pone.009361 7.g004 
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Figure 5. Arrival times, in ms, of the waves after point 
stimulation at the epicardial surface for a large anisotropy 
ratio Z?t:Z?2 = 1:0.1 1 1. The notation is the same as in Fig. 4. 
doi:1 0.1 371 /journal. pone.0093617.g005 

fact that the red lines are always located below the black lines 
shows that the increase of the fiber rotation angle results in faster 
wave propagation. 

Also, all solid lines correspond to Di\D2 — 1:0. 1 1 1, while dashed 
lines correspond to Di'.D2 — 1:0.25. If we now compare the solid 
and dashed lines of the same color (third column in Fig. 9), we see 
that the red lines are closer to each other than the black lines. This 
indicates that in presence of higher fiber rotation (the red lines) the 
decrease of D2 (i.e., solid vs. dashed) has a smaller effect on the 
arrival time. This once again illustrates that anisotropy is 
compensated for by the rotation of the fibers. 

In addition, we see in Fig. 9 that the red lines always have a less 
steep slope than the corresponding black lines. As going from black 
to red shows an increase in fiber rotation, we can conclude that the 
increase in rotation makes propagation faster in all cases. 

Scroll wave dynamics 

We have also studied scroll wave dynamics for the same values 
of anisotropy and fiber rotation. We generated a single scroll wave 
located approximately at the middle between the apex and the 
base of the ventricle and studied its behavior for different model 
parameters yo? yi and Di'.D2. We found that anisotropy 
substantially affects the dynamics of scroll waves. In all cases, 
the increase in the fiber rotation angle results in a decrease in the 
period of scroll wave rotation (Fig. 10). We see that for 
Di'.D2— 1:0.25, when the fiber rotation angle is increased from 
16° to 174°, the period drops significantly, from 277 to 257 ms. 
For Di'.D2— 1:0.111, we see similar dependency. In addition, the 
period for the same rotation angle for Di'.D2 — 1:0. 1 1 1 was slightly 
longer than for Di'.D2 = 1:0.25. 



Di:D2=l:0.111 



16° 



69° 



133° 



174° 



1163 



Figure 6. Arrival times, in ms, of the waves after point 
stimulation at the endocardial surface for a large anisotropy 
ratio Z?i:Z?2 = 1:0.1 1 1. The notation is the same as in Fig. 4. 
doi:10.1371/journal.pone.0093617.g006 

The scroll wave dynamics was also substantially affected by the 
anisotropy. In all cases, we observed a drift of the filament (Fig. 11). 
The drift always had two components, both in the vertical (i/^) and 
circumferential {(j)) directions. The total velocity of drift (Fig. 12 A) 
was very small, about 1 mm/s, which is about 0.2 mm per 
rotation; however, the drift was monotonic and persistent. The 
value of velocity had no clear relationship with the rotation angle; 
for Di'.D2 — 1:0.25, we see some tendency for velocity decrease 
with an increase in fiber rotation, while for Di'.D2— 1:0.111, the 
dependency is strongly non-monotonic and it is maximal for the 
intermediate values of the fiber rotation angle. The drift direction 
can be seen from the sign of the vertical and horizontal 
components of the velocity (Figs. 12B, C). Here again, the 
direction is affected by the rotation angle; however, we also did not 
find any clear tendency for either drift to the apex or base of the 
heart depending on the rotation angle. 

For D1.D2 — 1:0.25, the initial scroll wave was always stable and 
did not break down to multiple scroll waves. For Di'.D2 — 1:0.1 1 1, 
we did observe formation of the additional sources of excitation. 
However, in most cases, they appeared simultaneously at a 
substantial distance from the initial filament and not as a result of 
filament buckling and breakup due to rotational anisotropy in the 
way it was reported in [33]. The onset of new sources had a clear 
correlation with the fiber rotation angle (Fig. 13). We did not 
observe any instabilities for small and big a (models 1 and 4 in 
Table 1, Figs. 13A, D), however, for intermediate and large values 
of a (Figs. 13B, C), we observed new sources, and their number 
increased with the increase of the fiber rotation angle (compare 
panels B and G). Note that cases presented in panels B and G 
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Figure 7. Arrival times, in ms, of the waves after point 
stimulation at the epicardial surface for an intermediate 
anisotropy ratio D-yiD^. = 1 :0.25. The notation is the same as in Fig. 4. 
doi:1 0.1 371 /journal. pone.0093617.g007 

correspond to the largest drift velocities of a scroll wave; thus, the 
onset of secondary filaments may be related to the filament drift. 

We have also studied change in filament shape over time. For 
small values of total fiber angle a, the filament remained 
transmural, nearly straight, and stable. This case is shown in 
Fig. 13 A. For intermediate a equal to 69° and 133°, the filament 
not only drifted faster but also deformed to a transmural S or L- 
shape (Fig. 13B). For larger values of a, the filament again had a 
nearly straight shape (Fig. 13D). 

Discussion 

In this paper, we used our recent anatomical model of the LV of 
the human heart using a special CS (y, ij/, (j)) which gives an explicit 
analytical map from a rectangular domain to the heart shape and 
fiber orientation field. This allowed us to represent the heart's 
geometry on a rectangular grid and explicitly write expressions for 
boundary conditions. This approach may be helpful for studies of 
any phenomena in which boundary effects are of great impor- 
tance. 

One important feature of our model is the possibility to change 
the properties of anisotropy. The most significant characteristic of 
LV anisotropy is the rotation of myocardial fibers through the 
myocardial wall. As shown in [45], relatively simple rule-based 
global models of LV myofiber directions yield no worse results 
than complicated image-based locally optimized models. Since any 
locally optimized model needs to be regular enough and the 
regularity requires smoothing and interpolation, no image-based 
model can be an untouched copy of a real heart. 

We can change the degree of the fiber rotation in a consistent 
way and study its effect on normal and abnormal wave 
propagation in the heart. In this paper, we investigated two main 
features of wave propagation using the detailed human ventricular 
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Figure 8. Arrival times, in ms, of the waves after point 
stimulation at the endocardial surface for an intermediate 
anisotropy ratio Z?i:Z?2= 1:0.25. The notation is the same as in Fig. 4. 
doi:10.1371/journal.pone.0093617.g008 

electrophysiological model: the distribution of effective excitation 
speed and the dynamics of transmural scroll waves. 

In our study of the effect of the fiber rotation and anisotropy on 
wave propagation, the initial stimulation area was located at the 
apex and on the lateral epi- and endocardium of the LV. We 
found that the rotation of myocardial fibers accelerates the spread 
of excitation waves in the heart, which was explicitly demonstrated 
using models with different fiber rotation angles. This acceleration 
of wave propagation was discussed in [32]; it occurs because the 
wave can propagate with maximal speed in more directions with a 
larger rotation angle, which results in an overall faster wave 
propagation. Note that if the rotation angle is 2n or more, 
propagation with maximal speed will be possible in any direction, 
and in the limit of a large medium, the arrival time will be the 
same as in an isotropic medium with the velocity determined by 
the velocity along the fiber [32]. We were able to demonstrate this 
in an anatomical setup in which we explicitly changed the rotation 
of the fibers, while in [32] such an assessment was made for a 
single anisotropy configuration. While Young and Panfilov 
represented the tissue as a simple rectangular 3-D slab with 
plane-parallel fibers [32], we adopt a fully 3-D fiber architecture 
together with an LV shape. In this paper, we particularly 
considered a more realistic ventricular architecture and morphol- 
ogy that includes the following features: 

1 . A more realistic method for the fiber rotation angle around the 
axes [27], so called "Japanese-fan arrangement" [29]; and 

2. A realistic change of the fiber rotation angle values due to the 
displacement of the transmural axis from the LV apex to the 
base [27]. 
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Figure 9. Arrival times, in ms, as a function of the distance from the stimulation point for the apical M), epicardial {B\, and 
endocardial (O stimulation. The distance on the horizontal axis is measured in ms as the arrival time of the wave in the isotropic model (see text 
for more details). The red lines represent numerical experiments for total rotation angle a = 174°; the black lines represent for a= 16°; and the blue 
lines represent isotropy. The solid lines correspond to the case Di:D2 = 1:0.1 1 1, while the dashed lines correspond to the case Di:D2 = 1:0.25. The 
vertical segments display minimal and maximal arrival times in each group of nodes. The average, min, and max arrival times are displayed on the 
leftmost panels for Di:D2 = 1:0.1 1 1 and in the middle column for Di:D2 = 1:0.25. The right column compares the average arrival times. 
doi:1 0.1 371 /journal. pone.0093617.g009 



As in [32], our model also shows the faster excitation 
propagation with an increase in the fiber rotation angle and a 
decrease in the anisotropy. 

A second set of results in this paper concerns the dynamics of 
transmural scroll waves. The negative correlation found for the 



rotation period versus fiber rotation angle in Fig. 10 is different 
from the observations in [46] made by Qu et al, who observed an 
increasing period with a faster fiber rotation rate. Note, however, 
that their simulations used in-plane fiber rotation, while we work 
with a 3-D ventricular geometry. Both complementary cases can 
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Figure 10. Scroll wave rotation period T, ms, as a function of total fiber rotation angle a, deg. 

doi:10.1371/journal.pone.0093617.g010 



be qualitatively explained on geometrical grounds. In [46], the 
period was only affected by fiber rotation, which causes negative 
intrinsic curvature 7Z of the associated curved space [32]. By its 
definition, negative geometrical curvature represents a saddle-like 
space with positive angular deficit. More specifically, in a space 
with curvature 7Z, the circumference C of a circle (ball) with small 
radius r amounts to [47] 



C^lnrl 1 



(10) 



As the spiral tip in each cross-section needs to travel along a larger 
closed path of length C before completing a period, negative TZ is 
expected to increase the rotation period. Therefore, the trend 
found in [46] can be expected if the fiber rotation is the main 
determinant of the rotation period. In our present model, 
however, a transmural filament consists of spiral waves' tips in 
different layers of constant depth y. These y-surfaces are sphere- 
like and therefore have positive TZ. Thus, under normal 
excitability, the sphericity of the LV will decrease the rotation 
period [48,49]. We conclude that in general, the rotation period of 



scroll waves in the heart may decrease or increase with increasing 
fiber rotation angle a, depending on the relative strengths of the 
competing effects of the fiber rotation rate and the extrinsic 
curvature (sphericity) of the LV cavity. 

Regarding the non-monotonic dependence of filament drift 
velocity versus total rotation angle, we first note that the creation 
of secondary filaments in the regime of intermediate a is consistent 
with the clinical or experimental picture of a "mother rotor" 
during cardiac arrhythmias [50,51]. In such a scenario, the 
primary filament (mother rotor) remains stable and creates 
secondary sources that further disturb the electrical excitation of 
the heart, leading to cardiac arrest. In our simulations, we also see 
a similar situation with a stable mother rotor, which was always 
sustained until the end of the simulation time, and secondary 
excitation sources induced at some distance from it. 

In addition to the direct simulation results detailed and 
discussed above, our anatomical model [27] coupled to detailed 
electrophysiological model [11] may prove useful in future studies 
for the following reasons. First, our model can be used to 
investigate the possible contribution of the LV geometry to the 
propagation of the excitation waves. In particular, in certain heart 
diseases (e.g. dilated and hypertrophic cardiomyopathy, eccentric 




r 

- 0 
-25 
-50 

-75 

-86 




40 



20 



Figure 11. Potential, mV, on the LV surface during scroll wave rotation (left) and tip trajectory for Z?i:Z?2= 1:0.1 1 1 (red line) and for 
Z?i:Z?2= 1:0.25 (black line) (right). The results are shown for model 2 (yo = 0.2, yi = 0.7, see text and Table 1 for details). 
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Figure 12. Velocity of scroll wave filament drift for the simulation of 8 s. Average filament velocity i/c, mm/s (A). Velocity components 
multiplied by 1000, per second: latitudinal component (B) and longitudinal component (C). 
doi:1 0.1 371 /journal.pone.009361 7.g01 2 



and concentric cardiac hypertrophy, etc., see chapter 8 in [52]), 
the shape of the ventricle becomes more spherical and the 
thickness of the wall also increases. Such changes in geometry can 
easily be accommodated in our model. In general, the study of the 
effects of the LV geometry on excitation seems to be of great 



importance because many cardiac pathologies tightly correlate 
with changes in the LV geometrical characteristics. The LV 
becomes more dilated near the apex and thicker near the base 
during stress-induced ("Takotsubo") cardiomyopathy, or transient 
apical ballooning syndrome (see chapter 8 in [52]). Such 
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Figure 13. Scroll wave filaments in the LV model. The anisotropy ratio is Di:D2 = 1:0.1 11. Panels A, B, C, D: models 1, 2, 3, 4 (see Table 1), fiber 
rotation angle in the LV wall increases from panel A to panel D. The epicarduim (semitransparent colored surface; color denotes height from the red 
base to the purple apex), the endocardium (opaque white meshy surface), and filaments (black lines and dots). 
doi:1 0.1 371 /journal.pone.009361 7.g01 3 



remodeling of the LV geometry might be mimicked in similar 
mathematical models via the fitting of the geometric parameters 
values to account for the LV shape of a particular pathology. 

The results of our simulations can be verified by direct 
measurements of wave propagation on the whole heart prepara- 
tions, such as [53]. Note, however, that another factor important 
for overall excitation of the heart is the Purkinje conduction 
system. In order to compare experiments with our simulations, 
such measurements should be performed after chemical ablation 
of the Purkinje system [54]. 

Another potentially interesting application of this approach may 
be its application in studies on the defibrillation of cardiac tissue, in 
which case tissue texture and boundary conditions are also of key 
importance [55] . However, a bi-domain representation of cardiac 
tissue must be used for defibrillation [55-57], which does not fall 
under our present scope. Nonetheless, the extension of our 
approach for such cases is straightforward. The formulae for the 
diffusion term (7) and for the boundary conditions (9) can be 
directly used for representation of the bi-domain equations in the 
special CS. Then, the finite difference problem can be formulated 
in the same way as in our case and can be solved using any existing 
method (see [58]). 

In this paper, we studied wave propagation due to point 
stimulation and scroll waves. Other important wave propagation 
regimes include various types of scroll waves and turbulent patterns 
[19,59-61]. It was shown that heterogeneity [62,63] and anisotropy 
[31,64,65] of the tissue are significant factors determining the 
dynamics of these sources. The effect of heterogeneity on the 
dynamics of spiral waves was also studied in a series of papers by 
Shajahan, Sinha and co-authors [30,41,55,67]. In particular, in [66] 



they showed that some changes in the position, size, and shape of a 
conduction inhomogeneity can transform a single rotating spiral to 
spiral breakup or vice versa. Since our model provides tools for 
changing anisotropic properties and allows one to add heterogene- 
ity, these effects can also be studied using our approach. 

Appendices 

A Construction of spiral surfaces 

We model myocardial sheets as surfaces filling the LV. The 
filling was obtained by rotation of one surface around the vertical 
LV symmetry axis. We call these surfaces "spiral" (see Fig. 2). 

We introduce the special CS (y, ij/, 0), which is linked with the 
cylindrical CS (see (1) and (2)). In this CS, the equation of the 
spiral surface is 



(/>(y,lA) = (^0+</'maxy. 



(A.l) 



where different values of (/>o allow us to obtain different spiral 
surfaces and (/>max is a constant of the model. 

The parametrical equation of a spiral surface in cylindrical CS is 

= (r/, + (l-y)/)(ecosiA + (l-e)(l- sirnA)), 
ziy,^) = (z, + (1 -y)h){l - sin lA) + (1 -y)h. 
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B The myofibers' equations 

The equations are (see Fig. 3): 

where different values of the parameter ^ G (0, 1) correspond to 
different myofibers, 

O G [arcsin Y,n— arcsin Y], 

^ess{.Pi^) is the exphcit equation of a spiral surface in the cylindrical 
CS. 

At each point (y, i/^, 0), 0<i/^<7r/2, a fiber orientation vector 
v = (x',j', z') is given by the formulae: 



, dp - 

x = -—cosd — psmd——, 



dp . . . dS^ 



(B.l) 



(B.2) 



C Laplacian in implicit curvilinear coordinates 

The Laplacian is an important term in the reaction-diffusion 
equation as it is responsible for the modeling of electrical wave 
spreading. It can be written as div(D grad Jj where / is the 
transmembrane potential and D is an anisotropic local diffusion 
matrix. 

Below, we calculate the Laplacian for anisotropic diffusion in 
the Cartesian and the special CS. 

C.l The Cartesian CS. For an arbitrary diffusion matrix 



IJ J ^ IJ ^ J 

In consideration of Eq. (6), one can write: 



div(D grad/) = (A -^2)-E%^-|^ 



dxj dxi 



dxi dxj 



C.2 The curvilinear CS. Here we deduce from Eq. (C.2) the 
proper form of the Laplacian in the curvilinear coordinates ^q, ^i, 
^2- Let us consider and / as functions of ^q, ^i, ^2 calculate 
three types of derivatives. First, one has 



, dz dp dz d^ 

^ ^Jp'm^^'m' 

where p is determined by Eq. (1), 

a) = 7r(l-y), 



(B.3) 



dvj 



dxj 



where 



dvi 
dxj 



dxj dxj 



^dik'dxj' 



(C.3) 



(C.4) 



^ = ((rb + yl) cot (yn) - ^ 



d(D K ' 



Secondly, we evaluate 



and thirdly. 



8xi ^8^j 8xi' 



(C.5) 



dxi dxj 



y d'f d^k 



(C.6) 



where 



dz _ Zfe + (1 —y)h 
d~p~ rt + {\-y)l 



r. 



1L_ 
dxi d(k 



54 di,j ' dxi ' 



{C.l) 



dz W, . , ^ reJ{zt + {\-y)h) \ 



where (^max is a dimensionless parameter affecting fiber twist. 



The difficulty now lies in the fact that the functions Xj{^^ define 
the ^f. only implicitly. To evaluate the necessary derivatives, we 
need the following matrices: 



i={Jij)= 



(C.l 
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dvi 



dvi 



\k_i TTk\ 



d'xk 



(C.9) 



(C.IO) 



The matrices are linked between themselves with the following 
relations: 



D.1,2 The epicardium. On the epicardium, Up = {Zb'^h) cos ij/, 
nz= — {rij-\- l)es, so one can write (D.2) as 



^ izb + A) COS lA - irb + l)e. = 0. (D.4) 
op oz 



Everywhere on the epicardium, except the apex, 
df (r* + /)•£, df 



dp (ziy-\-hy COS il/ dz' 



(D.5) 



W = SJ, 



(C 11) apex, cos i/^ = 0, es = €>0, so the boundary condition 

looks like ^ ~ ^• 

D.1.3 The endocardium. On the endocardium, np = Zb cos 
ij/, nz= —rtjEs, so (D.2) can be written as 



We substitute (C.3), (C.5), (C.7) to (C.2) and get: 



div(D grad/) = Y,Pk' JT + Y.^^^' 



where 



Pk = D2 tr - 02)- ((Jv)^-tr(S J)+(JSJv)^+vTt^v) , (C. 12) 



qi^i are elements of matrix Q; 



Q = JDJ^. 



(C.13) 



— -Zb cosy/-— -ncs = 0. 
dp dz 



Everywhere on the endocardium, except the apex, 
df nes df 



dp Zb cos \j/ dz ' 



(D.6) 



(D.7) 



On the apex of the endocardium, like the epicardium, the 

boundary condition is =0. 

dz ^ 
D.2 Isotropic case, special CS. In the special CS, we use 



df 

and ^r-. Let us take into account that 
dy/ 



dy 



D Boundary conditions 

Let n be a normal vector to one of the LV boundary surfaces. 
For outer domain boundaries, we use a no-flux condition on the 
current. 

D.l Isotropic case, cylindrical CS. One can write the 
boundary condition as 



dn 



(p,0,z) = O. 



(D.l) 



By the definition of directional derivative and since the normal 
vector to the LV boundary lies in our problem always in the 
corresponding meridional half-plane. 



(D.2) 



where Up and are the normal vector components in the 
meridional half-plane. 

D.1.1 The equator. On the equator, np = 0, so (D.2) reduces 

to 



dz 



=0. 



(D.3) 



and 



dp dy dp dij/ dp 



dz dy dz dij/ dz 



(D.8) 



(D.9) 



D.2.1 The equator. Formula (D.3) can be rewritten as 
df (r, + (l-y)/)(l-6) df 



dij/ 



I 



dy 



(D.IO) 



D.2. 2 The epicardium. Formula (D.4) becomes 

df l{rb-\-l)ec^s — {zb-\-h)hsm\jj cos \jj df 
dy (zb + hf cos2 + (rz, + ife^ d\j/ ' 

D.2. 3 The endocardium. Formula (D.6) yields 

df Irb^cCs — ^bh sin i/^ cos ^ df 
dy z\ cos^ y\f + r^ej d\\i ' 



(D.ll) 



(D.12) 
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Note that if e > 0, the two last formulae cannot have zero values 
in the denominators and can be directly applied to the LV apex 
points. 

D.3 Anisotropic case, special CS. The boundary condition 



nD grad/ = 0, 



(D.13) 



with n, the normal to the LV surface. 

D.3.1 The epi- and endocardium. Let us write (D.13) in 
detail: 



ij \ k 



d^u Sxj 



= 0; 



we substitute (C.5) to this equation: 



= 0. 



Let us express ^r-: 
dy/ 



3(1) 



We can calculate the derivatives of y, 0, xj/ by x, j, z'- 



COS 6 sin 6 

yx=yppx = — Y~ ' yy=yppy = — Y~ ' 



yz = - 



sin(^ COS0 



n^D 1^ + n^D ^ + n^D ^ = 0. (D. 14) 

cy dy/ d(j) 

Here, the J^'*^'*^ are columns of derivatives of these special variables 
by Cartesian coordinates (see (C.8)): 



J^= dy/dy , 
\dyldz) 

and so on. 

D.3. 2 The equator. Let us write (D.13) detailed as the 
following: 

n^D grad/), +«^(Z) gradA grad/), = 0. (D. 15) 



Vector n is coUinear to the Oz axis, so /z^^.^ /Zj; = 0. The equation 
(D.15) goes over 

(Z)grad/), = 0. 



Writing down the matrix product: 



here, = r^+(l -y)/, Zm = Zb+i}-y)h- 
So 



¥ _ 



(D.16) 



(- [D^^x + D^^y)z^+D^^p^{\ -£)) ^ +z^/-(-i)20 sin<^+7)2i cos<^) ^ 
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